Name, Vorname:
In diesem Versuch
3. Filterexport und -simulation
Abbildungen in diesem Notebook wurden konvertiert mit https://www.base64-image.de/ und in den HTML-Code eingebettet.
Für diesen Versuch benötigen Sie die Software pyfda (Python Filter Design and Analysis), die Sie per pip install pyfda installieren können (wenn Sie eine Python-Installation auf Ihrem Rechner haben), alternativ gibt es unter https://github.com/chipmuenk/pyfda Binaries für Windows und für Ubuntu 20.10 (funktioniert vermutlich auch mit anderen halbwegs aktuellen Distros).
Nach dem Praktikumsversuch exportieren Sie das Notebook mit Textantworten, Codezellen und Plots als HTML (File -> Export Notebook As ... -> Export Notebook to HTML) und reichen es in Moodle ein.</div>
import os, sys
module_path = os.path.abspath(os.path.join('..')) # append directory one level up to import path
if module_path not in sys.path: # ... if it hasn't been appended already
sys.path.append(module_path)
import dsp_fpga_lib as dsp
dsp.versions() # print versions
%matplotlib inline
import matplotlib.pyplot as plt
size = {"figsize":(12,5)} # Plotgröße in Inch
import numpy as np
import scipy.signal as sig
import wave
from IPython.display import Audio, display
#-----------------------------------------------------------------------------
def wav2np(filename):
""" Read the wav-file and convert it to a one or two-dimensional numpy array,
depending on the number of channels.
Properties of the WAV-file are stored as function attributes (evil)
"""
wf = wave.open(filename,'rb')
wav2np.N_CH = wf.getnchannels() # number of channels
wav2np.W = wf.getsampwidth() # wordlength per sample in bytes
wav2np.N_FR = wf.getnframes() # number of frames
wav2np.f_S = wf.getframerate() # sample (frame) rate
print("{0} channels with {1} frames of {2} bytes and f_S = {3} Hz.".format(wav2np.N_CH, wav2np.N_FR, wav2np.W, wav2np.f_S))
if wav2np.W == 2:
samples_in = np.frombuffer(wf.readframes(-1), dtype=np.int16) # read wav data as 16 bit integers, R and L samples are interleaved
elif wav2np.W == 1:
samples_in = np.frombuffer(wf.readframes(-1), dtype=np.int8) # read wav data as 8 bit integers, R and L samples are interleaved
else:
raise TypeError("Unknown data format: {0} bytes".format(wav2np.W))
samples = np.array([samples_in[idx::wav2np.N_CH] for idx in range(wav2np.N_CH)], dtype=np.int32) # deinterleave channels to numpy array N_CHAN x N_FRAMES
return samples
Python version: 3.10.9 Numpy: 1.23.5 Scipy: 1.10.0 Matplotlib: 3.7.0 module://matplotlib_inline.backend_inline

Auf der linken Seite wählen und entwerfen Sie Filter unter den Tabs "Specs" (Filterentwurf), "b,a" (Eingabe von Filterkoeffizienten), "P/Z" (Eingabe von Polen und Nullstellen). Auf der rechten Seite schauen Sie sich Betrags- und Phasengang, Gruppenlaufzeit, P/N Diagramm und die transiente Antwort auf verschiedene Stimuli an.
Frequenzen können in unterschiedlicher Form eingeben und angezeigt werden:
Amplituden werden in V, W oder dB eingegeben und angezeigt.
Filterentwurfsalgorithmen benötigen unterschiedliche Eingabeparameter. Wenn "Order Minimum" angeklickt ist, versucht ein Algorithmus die minimale Ordnung zu den gegebenen Filterspezifikationen zu berechnen, halbwegs genau funktioniert das aber nur bei Equiripple (FIR) und den meisten IIR-Filtern. Bei pyfda können Sie nach dem Entwurf mit minimaler Ordung das Häkchen entfernen und dann selber Entwurfsparameter bearbeiten.
Nicht benötigte Parameter sind entweder deaktiviert oder ausgeblendet. Manche Parameter sind ausgegraut, aber bearbeitbar: Diese Parameter beeinflussen nicht den Filterentwurf, sie dienen aber z.B. dazu um das Toleranzschema anzuzeigen.
Ein Audiotrack von einer CD ($f_S = 44,1 \,\text{kHz}$) soll tiefpassgefiltert werden. Die Filterspezifikationen sind:
Es sollen die folgenden Filtertypen getestet werden:
Ermitteln Sie die folgenden Filtereigenschaften und fassen Sie sie in einer Tabelle zusammen mit den Spalten:
| Filtertyp | IIR/FIR? | Ordnung N | Gruppenlaufzeit τg | Mult. L pro Sample | Max. Samplerate [GHz] | Vorteil | Nachteile |
|---|---|---|---|---|---|---|---|
| Fourier-Approximation / Least-Square(Windowed FIR mit Boxcar-Fenster) |
FIR | 164 | 1.84 ps | 164 | 1.83 | geringe Nichtlinearität | hohe Rippel / zu geringe Steilheit |
| Windowed FIR mit Kaiser-Fenster | FIR | 268 | 30.2 ns | 268 | 1.12 | keine Rippel im Sperrbereich / hohe Däpfung im Sperrbereich | mittlerer Rippel im Sperrbereich |
| Equiripple | FIR | 163 | 18.4 ps | 163 | 1.8 | geringe Nichtlinearität | hoher Rippel |
| Butterworth | IIR | PyFDA -> max. 20 | 1.46 ms | 23 | 13.0 | kein Rippel / geringe Nichttlinearität | geringe Steilheit |
| Chebychev 1 | IIR | 14 | 4.75 ms | 17 | 17.6 | mittlere Steilheit | hoher Rippel im Durchlassband / hohe Nichttlinearität |
| Chebychev 2 | IIR | 14 | 1.48 ms | 17 | 17.6 | mittlere Steilheit | hoher Rippel im Sperrbereich / hohe Nichttlinearität |
| Elliptisch (Cauer) | IIR | 7 | 2.31 ms | 10 | 30 | geringe Ordung | hoher Rippel |
Welche der Filter sind FIR, welche IIR?
Vergleichen Sie die benötigte Ordnung N. Für nicht alle der Filter gibt es gut funktionierende Designalgorithmen, Sie müssen dann von Hand die Filterordnung und / oder Eckfrequenzen so lange anpassen, bis alle Spezifikationen erfüllt sind. Im Tab "Info" sehen Sie, ob Spezifikationen verletzt werden.
Welche Ordnung benötigen Sie, wenn Sie die Abtastrate vor der Filterung auf die Hälfte reduzieren? (Tipp: Aktivieren Sie das "Schloss" neben der Samplingfrequenz, bevor Sie f_S ändern.) Es genügt, wenn Sie das für ein Filter ausprobieren. Was müssen Sie beachten wenn Sie die Samplerate reduzieren wollen?
Welche Gruppenlaufzeit τg haben die Filter bei $f = 2 \,\text{kHz}$? Schauen Sie sich im y[n] an, wie die Signalform eines 500 Hz Sägezahnsignals durch das Filter verformt wird.
Wieviele Multiplikationen L pro Sample werden benötigt? Berücksichtigen Sie dabei die Filterstruktur (Optimierung durch Linearphasigkeit möglich?), IIR Filter sollen in einer kaskadierten SOS-Struktur realisiert werden - siehe folgende Abbildungen.
Welche maximale Samplerate kann man auf einem DSP mit max. $300 \cdot 10^6$ Multiplikationen / s verarbeiten? Zusätzliche Informationen: How fast are DSPs?
Fügen Sie den Plot des Betragsgangs mit Phasengang (Checkbox "Phase") ein.
| Direktform | Lin. phasige Form | Kaskadierte Second Order Sections |
Windowed FIR mit Kaiser-Fenster:
Equiripple:
Butterworth:
Chebychev 1:
Chebychev 2:
Elliptisch:

Im folgenden sollen Sie Audiofiles mit dem Equiripple und dem elliptischen Filter filtern. Dazu können Sie aus pyfda Filter in verschiedenen Formaten exportieren:
a = [...]) in die nächste Codezelle einfügen.Sie können Daten im Zeitbereich filtern mit
y = sig.lfilter(b,a,x)
b und a sind die Koeffizientenvektoren (nicht-rekursiver / rekursiver Teil), x ist das Eingangssignal. y hat die gleichen Dimensionen wie x (es können auch Stereofiles verarbeitet werden). Was müssen Sie bei einem FIR-Filter für a angeben?
Wenn die Filterbeschreibung als Second-Order Sections vorliegt, kann auch der folgende Befehl benutzt werden (allerdings unterstützt pyfda noch keinen Export im sos - Format):
y = sig.sosfilt(sos,x)
Anhören können Sie sich die gefilterten Audiofiles wieder mit display(Audio(data=y, rate=wav2np.f_S)).
Sie können auch eigene WAV-Files in Ihr medien Verzeichnis hochladen und testen, bei großen Audiofiles sollten Sie zunächst nur ein paar Sekunden testhören (mit x[:N_samples] bzw. x[:, :N_samples]) bei Stereo-Samples. Was bewirkt die seltsame Syntax bei Stereosignalen?
Hören Sie einen Unterschied zwischen dem Orginal und den beiden unterschiedlichen Filtern? Woran liegt das?
x_w = wav2np("../medien/87778__marcgascon7__vocals.wav")
#Equiripple
b_FIR = [0.0006123033469102477,0.0005231312265310316,0.000561083771645731,0.00041518098768839355,2.827711686114232e-05,-0.0006169795177762187,-0.0014804911786779923,-0.002458026960357289,-0.0033955247678761798,-0.004109019417673139,-0.004428667554745472,-0.004237079147342417,-0.0035082858849691473,-0.002325963428206755,-0.0008793225682988921,0.0005669966191960529,0.001727314146414968,0.002358436650440632,0.0023225886047265588,0.0016281694762668173,0.00044002268486894485,-0.0009495414647988306,-0.00218391145568962,-0.002926600347304991,-0.002948523144002793,-0.0021936643092102604,-0.0008066063399116822,0.0008911226757635836,0.002471991422020751,0.0035073292743326347,0.0036806624276502827,0.0028800941086848594,0.0012444616340492698,-0.0008540601244899048,-0.0028900585715064993,-0.0043130326782202496,-0.004692235633235503,-0.00384162591921165,-0.0018896284603115599,0.0007297360484630355,0.0033677431056925913,0.005317244868695226,0.005996942088874641,0.005117766406001423,0.0027847237841484847,-0.0004979172934868248,-0.0039276848931389715,-0.0065934641947688385,-0.0077099517214559606,-0.006839091317099379,-0.004037685130741728,0.00011364684868242633,0.004618634275870731,0.008293677454136077,0.010066529782336353,0.009273031150922878,0.005873052144054021,0.0005221554755543486,-0.005533741821921171,-0.010729128825649713,-0.013560660832545542,-0.012991437230981336,-0.008775651715744494,-0.00161273198833575,0.0069258868738092455,0.014701407248243398,0.01949170936897492,0.019552878251152542,0.014129162521417389,0.003783120765910111,-0.00955025960828738,-0.02282413158635936,-0.03242421406707713,-0.03489768645270763,-0.027719185939259615,-0.009927599299774674,0.017513404923713328,0.051725839878917296,0.08831767123202569,0.12210814306160755,0.14804243901826264,0.1621149051862032,0.1621149051862032,0.14804243901826264,0.12210814306160755,0.08831767123202569,0.051725839878917296,0.017513404923713328,-0.009927599299774674,-0.027719185939259615,-0.03489768645270763,-0.03242421406707713,-0.02282413158635936,-0.00955025960828738,0.003783120765910111,0.014129162521417389,0.019552878251152542,0.01949170936897492,0.014701407248243398,0.0069258868738092455,-0.00161273198833575,-0.008775651715744494,-0.012991437230981336,-0.013560660832545542,-0.010729128825649713,-0.005533741821921171,0.0005221554755543486,0.005873052144054021,0.009273031150922878,0.010066529782336353,0.008293677454136077,0.004618634275870731,0.00011364684868242633,-0.004037685130741728,-0.006839091317099379,-0.0077099517214559606,-0.0065934641947688385,-0.0039276848931389715,-0.0004979172934868248,0.0027847237841484847,0.005117766406001423,0.005996942088874641,0.005317244868695226,0.0033677431056925913,0.0007297360484630355,-0.0018896284603115599,-0.00384162591921165,-0.004692235633235503,-0.0043130326782202496,-0.0028900585715064993,-0.0008540601244899048,0.0012444616340492698,0.0028800941086848594,0.0036806624276502827,0.0035073292743326347,0.002471991422020751,0.0008911226757635836,-0.0008066063399116822,-0.0021936643092102604,-0.002948523144002793,-0.002926600347304991,-0.00218391145568962,-0.0009495414647988306,0.00044002268486894485,0.0016281694762668173,0.0023225886047265588,0.002358436650440632,0.001727314146414968,0.0005669966191960529,-0.0008793225682988921,-0.002325963428206755,-0.0035082858849691473,-0.004237079147342417,-0.004428667554745472,-0.004109019417673139,-0.0033955247678761798,-0.002458026960357289,-0.0014804911786779923,-0.0006169795177762187,2.827711686114232e-05,0.00041518098768839355,0.000561083771645731,0.0005231312265310316,0.0006123033469102477]
a_FIR = [1.]
#Elliptisch
b_IIR = [0.001710952004172366,-0.005755838859144107,0.008363025285148405,-0.004112746893907044,-0.004112746893907043,0.008363025285148403,-0.005755838859144107,0.001710952004172366]
a_IIR = [1.0,-6.10402966946488,16.40268768530511,-25.100822235115835,23.595177173520497,-13.615049072423911,4.4641836431418325,-0.6417367418902715]
#Signal Filtern
y_FIR = sig.lfilter(b_FIR, a_FIR, x_w)
y_IIR = sig.lfilter(b_IIR, a_IIR, x_w)
print("Orginal")
display(Audio(data=x_w, rate=wav2np.f_S))
print("-----------------------------------------------------------")
print("FIR: Equiripple")
print(y_FIR.dtype)
display(Audio(data=y_FIR, rate=wav2np.f_S))
print("-----------------------------------------------------------")
print("IIR: Elliptisch")
print(y_IIR.dtype)
display(Audio(data=y_IIR, rate=wav2np.f_S))
2 channels with 349952 frames of 2 bytes and f_S = 44100 Hz. Orginal
----------------------------------------------------------- FIR: Equiripple float64
----------------------------------------------------------- IIR: Elliptisch float64
#fig, ax = plt.subplots(**size)
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(35,8))
l = 4000
start = 20000
ax[0].set_title("FIR-Filter: Equiripple")
ax[0].plot(np.arange(l), x_w[0][start:start+l], label="Orginal")
ax[0].plot(np.arange(l), y_FIR[0][start:start+l], label="y_FIR")
ax[0].legend()
ax[1].set_title("IIR-Filter: Elliptisch")
ax[1].plot(np.arange(l), x_w[0][start:start+l], label="Orginal")
ax[1].plot(np.arange(l), y_IIR[0][start:start+l], label="y_IIR")
ax[1].legend()
<matplotlib.legend.Legend at 0x1b3ef61a380>
(c) 2016 - 2021 Prof. Dr. Christian Münker
This jupyter notebook is part of a collection of notebooks on various topics of Digital Signal Processing. The latest version can be found at https://github.com/chipmuenk/dsp.
This notebook is provided as Open Educational Resource. Feel free to use it for your own purposes. The text is licensed under Creative Commons Attribution 4.0, the code of the IPython examples under the MIT license. Please attribute the work as follows: Christian Münker, Digital Signal Processing - Vorlesungsunterlagen mit Simulationsbeispielen, 2020.